Spin waves in paramagnetic BCC iron: spin dynamics simulations 
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Large scale computer simulations are used to elucidate a longstanding controversy regarding the exis- 
tence, or otherwise, of spin waves in paramagnetic BCC iron. Spin dynamics simulations of the dynamic 
structure factor of a Heisenberg model of Fe with first principles interactions reveal that well defined peaks 
persist far above Curie temperature Tc. At large wave vectors these peaks can be ascribed to propagating 
spin waves, at small wave vectors the peaks correspond to over-damped spin waves. Paradoxically, spin 
wave excitations exist despite only limited magnetic short-range order at and above Tc. 
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For over three decades, the nature of magnetic excita- 
tions in ferromagnetic materials above the Curie tempera- 
ture Tc has been a matter of controversy amongst experi- 
mentalists and theorists alike. Early neutron scattering ex- 
periments on iron suggested that spin waves were renor- 
malized to zero at |1]; however, in 1975, using unpo- 
larized neutron scattering techniques, Lynn at Oak Ridge 
(ORNL) reported |2] that spin waves in iron persisted as 
excitations up to the highest temperature measured (lATc), 
and no further renormalization of the dispersion relation 
was observed above T^,. 

Experimentally, it was challenged primarily by Shirane 
and collaborators at Brookhaven (BNL) |3]. Using polar- 
ized neutrons, they reported that spin wave modes were 
not present above Tc and suggested that the ORNL group 
needed polarized neutrons to subtract the background scat- 
tering properly. Utilizing full polarization analysis tech- 
niques the ORNL group subsequently confirmed their ear- 
lier work and, in addition, they analyzed data from both 
groups and concluded that their resolution was more than 
an order of magnitude better than that employed by the 
BNL researchers [41. Moreover, angle -resolved photoe- 
mission studies jflQl suggested the existence of magnetic 
short-range order (SRO) in paramagnetic iron and that this 
could give rise to propagating modes. Theoretically, SRO 
of rather long length scales (25 A) was postulated to exist 
far above Tc |7, 8] and a more subtle kind was proposed 
later |9]. Contrarily, it was also suggested that above Tc, 
all thermal excitations are dissipative ifloLITTIl . To further 
complicate matters, analytical calculations for a Heisen- 
berg model of iron, with exchange interactions extending 
to fifth-nearest neighbors and a three pole approximation 
fl2ll . did not reproduce the line shape measured by either 
experimental group mentioned above. In addition, Shastry 
II3I1 performed spin dynamics (SD) simulations of a near- 
est neighbor Heisenberg model of paramagnetic iron with 
8192 spins and showed some plots of dynamic structure 
factor S{q,uj) with a shoulder at nonzero u for some q. 
It was explained to be due to statistical errors instead of 



propagating modes. 

With new algorithmic and computational capabilities, 
qualitatively more accurate SD simulations can now be per- 
formed. In particular, it can follow many more spins for 
much longer integration time. We use these techniques 
and a model designed specifically to emulate BCC iron and 
have been able to unequivocally identify propagating spin 
wave modes in the paramagnetic state, lending substantial 
support to Lynn's |2] experimental findings. Interestingly, 
spin waves are found despite only limited magnetic SRO. 

To describe the high temperature dynamics we use a 
classical Heisenberg model TC = —(1/2) X^^^^r' "^r.r'Sr • 
Sf / , for which the exchange interactions, Jr.rs are obtained 
from first principles electronic structure calculations. For 
Fe this is a reasonable approximation since the size of the 
magnetic moments associated with individual Fe-sites are 
only weakly dependent on the magnetic state lfl4ll and by 
including interactions up to fourth nearest neighbors it is 
possible to obtain a reasonably good T^. 

Large scale computer simulations using SD techniques 
to stiidv the dynamic properties of Heisenberg ferromag- 
nets II5I1 and antiferromagnets 1 16] have been quite effec- 
tive, and the direct comparison of RbMnFa SD simulations 
with experiments was especially satisfying |16]. We have 
adopted these techniques and used Lx Lx L BCC lattices 
with periodic boundary conditions and L = 32 and 40. At 
each lattice site, there is a three-dimensional classical spin 
of unit length (we absorb spin moments into the definition 
of the interaction parameters) and each spin has a total of 
50 interacting neighbors. We use interaction parameters, 
Ji, for the T = ferromagnetic state of BCC Fe calculated 
using the standard formulation llVll and the layer-KKR 
method Iil8ii1 . The calculated values are Ji = 36.3386 
meV, J2 = 20.6520 meV, J3 = -1.625962 meV, and 
J4 = -2.39650 meV. 

In our simulations, a hybrid Monte Carlo method was 
used to study the static properties and to generate equi- 
librium configurations as initial states for integrating the 
coupled equations of motion of SD At Tc and for 
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L = 32, the measured nonlinear relaxation time in the 
equilibrating process and the linear relaxation time be- 
tween equilibrated states for the total energy and for the 
magnetization |20] are both smaller than 500 hybrid steps 
per spin. We discarded 5000 hybrid steps (for equilibra- 
tion) and used every 5000*^* hybrid step's state as an ini- 
tial state for the SD simulations. For the J^'s used here, 
Tc = 919(l)i^, which is slightly smaller than the exper- 
imental value T^^^ = \QA2>K. The equilibrium magne- 
tization |m| = {l/N)\ Sri ~ (1 - T/T^y/^ in the 
vicinity of T^, and this is in agreement with experiments. 
The SD equations of motion are 
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r' Sr' is an effective field at site 



r due to its interacting neighbors. The integration of the 
equations determines the time dependence of each spin and 
was carried out using an algorithm based on second-order 
Suzuki-Trotter decompositions of exponential operators as 
described in jSlll . The algorithm views each spin as under- 
going Larmor precession around its effective field H^ff, 
which is itself changing with time. To deal with the fact 
that we are considering four shells of interacting neighbors, 
the BCC lattice is decomposed into sixteen sublattices. 
This algorithm allows time steps as large as 5t = 0.05 (in 
units of to = Ji^)- Typically, the integration was carried 
out to t„,ax = 20000(^t = lOOOto- 

The space- and time-displaced spin-spin correlation 
function C"''(r — r',t) and the related dynamical structure 
factor, S^Ja, w), are fundamental in the study of spin dy- 
namics ll23l and are defined as 
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where k = x, y or z and the angle brackets (• • • ) denote 
the ensemble average, and 
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where q and uj are momentum and energy (E oc to) trans- 
fer respectively. It is 5''^(q, to) that was probed in the neu- 
tron scattering experiments discussed earlier. 

By calculating partial spin sums 'on the fly' ifT^ . 
it is possible to calculate 5''^(q, w) without storing a 
huge amount of data associated with each spin config- 
uration. Because L is finite, only a finite set of q 
values are accessible: q = 27rn^/(La) with = 
±1, ±2, . . . , zizL for the {q, 0, 0) and (q, q, q) directions 
and riq = ±1, ±2, . . . , ibL/2 for the (q, q, 0) direction, 
(a is lattice constant.) For T > T^, the ensemble average 
in Eq. 121 was performed using at least 2000 starting config- 
urations. We average S''^(q, cj) over equivalent directions 
and this averaged structure factor is denoted as ^(q, uj). 

In Fig.^we show the frequency dependence of ^(q, cj) 
obtained for four different temperatures around T^. These, 
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FIG. 1: Calculated energy dependence of ^(q, oj) at q = 
7r/a(l,0,0) and for T = 0.95Tc (700 runs), l.QTc (2000 runs), 
l.lTc (2240 runs), and 1.2Tc (2240 runs) for L = 32. The solid 
lines are fits to the data as explained in text. Error bars are shown 
at a few typical points. 



so called, constant-q scans are for q = 7r/a(l,0,0) (|q| = 
1.09 A^^), which is half way to Brillouin Zone boundary. 
At 0.95rc, ^(q, u) already has a 3-peak structure: one 
weak central peak at zero energy and two symmetric spin 
wave peaks (we only show data for liJ > since the struc- 
ture factor is symmetric about u = 0). Note that the spin 
wave peaks are already quite wide. As T goes to T^, and 
above, the central peak becomes more pronounced. In ad- 
dition, the spin wave peaks shift to lower energies, broaden 
further and become less obvious, however they still persist. 
This 3-peak structure at high temperatures is in contrast to 
the 2-peak spin wave structure found at low temperatures. 
In the neutron scattering from ^'*Fe(12%Si) experiments 
|[4i], Mook and Lynn also noticed a central peak, but could 
not decide whether it was intrinsic to pure iron or a result 
of alloying of silicon. 

In general, constant-q scans are isotropic in the (g, 0, 0), 
(g, q, 0), and (q, q, q) directions. For very small |q|, there 
is only a central peak in the scans (as is expected) and the 
3-peak structure only develops for larger |q|. We fit the 3 
peaks in S'(q, to) using different fitting functions and found 
the best results with either a Gaussian central peak plus two 
Lorentzian peaks at zbojo: 

S{ci,uj) = G + L+ + L_, (4) 

or a Gaussian central peak plus two additional Gaussian 
peaks at iwo- 

S{q,u) = G + G+ + G., (5) 

where G = Icexp{—uj'^/ujl), L± = IqujI/{{uj =f uJoY + 
Lul),andG± = Ioea;p(— (cj=pu;o)^/wi)- Formoderate |q| 
the results are fit best with Eq.|3 while Eq.lsjworks better 
at larger |q|. In Fig.|2|we show, for T = T^, the results of 
fitting constant-q scans at |q| = 0.48 A^^ and |q| = 1.16 
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FIG. 2: Fits to 5'(q,cj) at T = Tc for two |q|-points along the 
((j,g,0) direction for L = 32. (a) |q| = 0.48 A'^ fit to Eq.El 
with 7c = 58.7, = 27.6 meV, Jq = 80.6, cjq = 21.4 meV, and 
wi = 11.3 meV; (b) |q| = 1.16 A^^ fit to Eq.|5lwith Ic = 2.49, 
= 193.3 meV, Iq = 3.02, luq = 175.3 meV, and lui = 61.2 
meV. The vertical scale in (b) is much smaller than that in (a). 
Error bars are shown at a few typical points. 



A"^ in the {q,q,0) ciirection. The |q| = 0.48 A"^ result 
fits well to Eq.|4]and has uji/uq < 1, i.e., the excitation 
lifetime is longer than its period and thus it can be regarded 
as a spin wave excitation. It should be noted that this |q| 
value is very close to that (0.47 A~^) for which Lynn found 
propagating modes in contradiction to the findings of the 
BNL group. At |q| = 1.16 A^^, the structure factor has 
much weaker intensity and fits best to Eq. |5] with a ratio 
u^i/loq that is even smaller than at |q| = 0.48 A~^. This is 
illustrative of the general conclusion that the propagating 
nature of the excitation modes is most pronounced at large 

|q|- 

Figure|3shows the dispersion relations obtained by plot- 
ting the peak positions, uJq, determined from the fits to 
^(q, uj) along the (q, q, 0) direction. Calculated disper- 
sion curves are shown at several temperatures in the fer- 
romagnetic and paramagnetic phases together with the ex- 
perimental results of Lynn |2]. To estimate errors, we fit- 
ted each constant-q scan several times by cutting off the 
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FIG. 3: Comparison of dispersion curves obtained in our sim- 
ulations (sim) with Lynn's experimental (exp) (Ref. results 
for the {q, q, 0) direction. Open symbols indicate excitations with 
mixed nature and are not due to spin waves (NSW). 



tail at slightly different Uj^ax to get an average cjq; these 
error bars are found to be no larger than symbols. In this 
figure, filled symbols indicate modes that are clearly propa- 
gating (uJi/uJo < 1) while open symbols indicate that, even 
though there are peaks at 7^ 0, the peaks have widths 
uJi > lOq. The calculated result for T = O.ST^ is very 
close to that from the experiments and propagating modes 
exist for very small |q|. For T > T^, our curves lie below 
the experiments's and soften with increasing temperatures, 
a property not seen in the experiments. One possibility de- 
serving of further study is that our use of temperature and 
configuration independent exchange interactions, in partic- 
ular those appropriate to the T = ferromagnetic state, 
breaks down at high temperatures when the spin moments 
are highly non-collinear. 

In our simulations we have equal access to constant-q 
scans and constant- ii^ scans; however, this is not the case 
in neutron scattering experiments. Because the dispersion 
curves of Fe are generally very steep, experimentalists usu- 
ally perform constant-E' scans. In Fig.|3]we show constant- 
E scans for several E values at T = LIT^ based on sim- 
ulations. Clearly, the constant-i? scans have two peaks 
(symmetric about |q| =0) that become smaller and wider 
and shift to higher |q| as increases. Peaks in constant-E' 
scans strongly suggest that SRO persists above |7]. 

The degree of magnetic SRO can be obtained directly 
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FIG. 4: T — l.lTc constant-_E scans along {q, q, 0) direction for 
E = 41.1 meV, 54.8 meV, 68.5 meV, and 96.0 meV with L = 40. 
Brillouin Zone boundary q^b = 155 A^^ in the direction. 

from the behavior of static correlation function C"^(r — 
r', 0) (i.e. Eq.Elwith t = 0), which can be calculated from 
the Monte Carlo configurations alone. For T = l.lTc we 
find a correlation length of approximately 2a (~ 6 neigh- 
bor shells), indicative of only limited SRO. Thus, in gen- 
eral, extensive SRO is not required to support spin waves. 
Moreover, inspection of Fig. |3lfor T = LIT^ shows that 
the point q > 0.77A^^, at which these peaks first cor- 
respond to propagating modes, is when their wavelength 
(A ~ 2a) first becomes the order of the static correlation 
length. 

In summary, our SD simulations clearly point to the ex- 
istence of spin waves in the paramagnetic state of BCC Fe 
and support the original conclusions of Lynn. Their sig- 
nature is seen as spin wave peaks in dynamical structure 
factor in constant-q and constant-£^ scans. Detailed analy- 
sis of the constant-q scans shows that the propagating na- 
ture of these excitations is clearest at large |q|, in agree- 
ment with experiment. This is also consistent with the re- 
quirement that their wavelength be the order of, or shorter 
than, the static correlation length. While the inclusion of 
four shells of first-principles-determined interactions into 
the Heisenberg model makes our results specifically relate 
to BCC Fe, we have also found spin waves in a Heisen- 
berg model containing only nearest neighbor interactions. 
In addition to elucidating the longstanding controversy re- 
garding the existence of spin waves above T^, these simula- 
tions also point to the important role that inelastic neutron 
scattering studies of the paramagnetic state can have in un- 
derstanding the nature of magnetic excitations, particularly 
when coupled with state-of-the-art SD simulations. 

We thank Shan-Ho Tsai, H. K. Lee, V. R Antropov, and 
K. Binder for informative discussions. Computations were 



performed at ORNL-CCS (www.ccs.ornl.gov) and NERSC 
(www.nersc.gov) using elements of the ^-Mag toolset 
which may be obtained at http://mri-fre.ornl.gov/psimag! 
Work supported by Computational Materials Science Net- 
work sponsored by DOE BES-DMSE (XT), BES-DMSE 
(GMS), NSF Grant No. DMR-0341874 (XT and DPL) 
and DARR\ (XT and TCS) under contract No. DE-AC05- 
00OR22725 with UT-Battelle LLC. 



* Electronic address: txp@uga.edul 

[1] M. F. Collins, V. J. Minkiewicz, R. Nathans, L. Passell, and 
G. Shirane, Phys. Rev. 179, 417 (1969). 

[2] J. W. Lynn, Phys. Rev. B 11, 2624 (1975). 

[3] J. P Wicksted, G. Shirane, and O. Steinsvoll, Phys. Rev. B 
29, R488 (1984); G. Shirane, O. Steinsvoll, Y. J. Uemura, 
and J. Wicksted, J. Appl. Phys. 55, 1887 (1984); J. R Wick- 
sted, P Boni, and G. Shirane, Phys. Rev. B 30, 3655 (1984). 

[4] H. A. Mook and J. W. Lynn, J. Appl. Phys. 57, 3006 (1985). 

[5] E. M. Haines, V. Heine, and A. Ziegler, J. Phys. F: Metal 
Phys. 15, 661 (1985); E. M. Haines, R. Clauberg, and 
R. Feder, Phys. Rev. Lett. 54, 932 (1985). 

[6] E. Kisker, R. Clauberg, and W. Gudat, Z. Phys. B 61, 453 
(1985). 

[7] V. Korenman, J. L. Murray, and R. E. Prange, Phys. Rev. B 

16, 4032 (1977); R. E. Prange and V. Korenman, Phys. Rev. 

B 19, 4691 (1978). 
[8] H. Capellmann, Solid State Commun. 30, 7 (1979); Z. Phys. 

B 35, 269 (1979). 
[9] V. Heine and R. Joynt, Europhys. Lett. 5, 81 (1988). 
[10] J. Hubbard, Phys. Rev. B 20, 4584 (1979); Phys. Rev. B 23, 

5974 (1981). 

[11] T. Moriya, J. Magn. Magn. Mat. 100, 261 (1991). 

[12] B. S. Shastry, D. M. Edwards, and A. P Young, J. Phys. C 

14, L665 (1981). 
[13] B. S. Shastry, Phys. Rev. Lett. 53, 1104 (1984). 
[14] A. J. Pindor, J. Staunton, G. M. Stocks, and H. Winter, J. 

Phys. F 13, 979 (1983). 
[15] K. Chen and D. P Landau, Phys. Rev. B 49, 3266 (1994). 
[16] S.-H. Tsai, A. Bunker, and D. P Landau, Phys. Rev. B 61, 

333 (2000). 

[17] A. I. Liechenstein, M. I. Katsnelson, V. P. Antropov, and 
V. A. Gubanov, J. Mag. Mag. Mater. 67, 65 (1987). 

[18] T. C. Schulthess and W. H. Butler, J. App. Phys. 83, 7225 
(1998). 

[19] P Peczak and D. P Landau, Phys. Rev. B 47, 14260 (1993). 

[20] D. P. Landau and K. Binder, A Guide to Monte Carlo Sim- 
ulations in Statistical Physics (Cambridge University Press, 
Cambridge, 2000). 

[21] M. Krech, A. Bunker, and D. P. Landau, Comput. Phys. 
Commun. Ill, 1 (1998). 

[22] S. W. Lovesey, Theory of neutron scattering from condensed 
matter (Clarendon, Oxford, 1984). 



